'''
 * @desc : svm支持向量机
 * @auth : TYF
 * @date : 2019/9/13 - 15:58

  在n维空间中能找到一个n-1维平面进行划分,找到分界面两边的样本点使得两点距离取值能够最大,则这些点称为支持向量,而n-1维分类超平面称为支持向量机

  n维空间平面,其中w(w1,w2,..)为法向量决定了超平面的方向,b为位移项决定了超平面与原点的距离
  w'x + b = 0

  点到平面距离公式
  r = |w'+b| / ||w||
  两个异类支持向量到分类平面的距离之和
  r = 2 / ||w||

  优化问题即找到一组w和b(即一个平面)使得异类支持向量间隔最大:
  设有分类平面
  f(x)= w'x + b
  找到一组w、b满足以下条件
  w'x(i) + b >= +1  y(i)=+1
  w'x(j) + b <= -1  y(j)=-1
  找到分类超平面的最大间隔,即最大化下式
  max (r)
  max ( 2 / ||w|| )

  拉格朗日乘子法（Lagrange Multiplier) 和KKT条件
  通常我们需要求解的最优化问题有如下几类：
  (i)无约束优化问题，可以写为:
      min f(x);  
      常常使用的方法就是Fermat定理，即使用求取f(x)的导数，然后令其为零，
      可以求得候选最优值，再在这些候选值中验证；如果是凸函数，可以保证是最优解。

  (ii)有等式约束的优化问题，可以写为:
      min f(x), 
      s.t. h_i(x) = 0; i =1, ..., n 
      对于第(ii)类的优化问题，常常使用的方法就是拉格朗日乘子法（Lagrange Multiplier) ，
      即把等式约束h_i(x)用一个系数与f(x)写为一个式子，称为拉格朗日函数，而系数称为拉格朗日乘子。
      通过拉格朗日函数对各个变量求导，令其为零，可以求得候选值集合，然后验证求得最优值。

  (iii) 有不等式约束的优化问题，可以写为：
      min f(x), 
      s.t. g_i(x) <= 0; i =1, ..., n
      h_j(x) = 0; j =1, ..., m
      对于第(iii)类的优化问题，常常使用的方法就是KKT条件。
      同样地，我们把所有的等式、不等式约束与f(x)写为一个式子，也叫拉格朗日函数，系数也称拉格朗日乘子，
      通过一些条件，可以求出最优值的必要条件，这个条件称为KKT条件。


  (a) 拉格朗日乘子法（Lagrange Multiplier)
  对于等式约束，我们可以通过一个拉格朗日系数a
  把等式约束和目标函数组合成为一个式子L(a, x) = f(x) + a*h(x),
  这里把a和h(x)视为向量形式，a是横向量，h(x)为列向量。
  然后求取最优值，可以通过对L(a,x)对各个参数求导取零，联立等式进行求取。
  (b) KKT条件
  对于含有不等式约束的优化问题，如何求取最优值呢？常用的方法是KKT条件，
  min f(x), 
 s.t. g_i(x) <= 0; i =1, ..., n
 h_j(x) = 0; j =1, ..., m
  同样地，把所有的不等式约束、等式约束和目标函数全部写为一个式子L(a, b, x)= f(x) + a*g(x)+b*h(x)，KKT条件是说最优值必须满足以下条件：
  1. L(a, b, x)对x求导为零；
  2. h(x) =0;
  3. a*g(x) = 0;
  求取这三个等式之后就能得到候选最优值。
  其中第三个式子非常有趣，因为g(x)<=0，如果要满足这个等式，必须a=0或者g(x)=0. 这是SVM的很多重要性质的来源，如支持向量的概念。


  序列最小优化(Sequential Minimal Optimization, SMO):
  上述二次规划问题规模正比于训练样本数，SMO算法工作原理是
  每次循环选择两个α进行优化处理，一旦找到一对合适的α就增大其中一个同时减小另外一个。
  选择的这两个α必须要在间隔边界之外。


'''

from __future__ import print_function
from numpy import *
import matplotlib.pyplot as plt


'''
样本示例
2.326297	0.265213	-1
3.634009	1.730537	-1
0.341367	-0.894998	-1
3.125951	0.293251	-1
2.123252	-0.783563	-1
0.887835	-2.797792	-1
7.139979	-2.329896	1
'''

'''
返回样本及其标签
'''
def loadDataSet(fileName):
    dataMat = []
    labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat, labelMat


'''
返回一个0~m之间，且不为i的整数
即已知第一个α为αi，随机找第二个αj，0<j<m
'''
def selectJrand(i, m):
    j = i
    while j == i:
        j = int(random.uniform(0, m))
    return j

'''
调整aj为H, L之间
'''
def clipAlpha(aj, H, L):
    if aj > H: #大于H则调整为H
        aj = H
    if L > aj: #小于L则调整为L
        aj = L
    return aj


"""
smoSimple
Args:
    dataMatIn    数据集
    classLabels  类别标签(+1或者-1)
    C   松弛变量(常量值0.6)，允许有些数据点可以处于分隔面的错误一侧。
        控制最大化间隔和保证大部分的函数间隔小于1.0这两个目标的权重。
           可以通过调节该参数达到不同的结果。
    toler   容错率（常量值0.001，是指在某个体系中能减小一些因素或选择对某个系统产生不稳定的概率。）
    maxIter 退出前最大的循环次数
Returns:
    b       模型的常量值
    alphas  拉格朗日乘子
"""
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    # 样本矩阵
    dataMatrix = mat(dataMatIn)
    # 标签转置成列
    labelMat = mat(classLabels).transpose()
    m, n = shape(dataMatrix)
    # 初始化b
    b = 0
    # 初始化α为全1
    alphas = mat(zeros((m, 1)))
    # 没有任何alpha改变的情况下遍历数据的次数
    iter = 0
    # 外循环用于确定最佳alphas对，每次确定好最佳alphas对后优化至所有alphas满足KKT条件(即收敛)更快
    # 这里简化版SMO算法忽略该过程，用迭代次数限制来代替，所以简化版SMO会比完整版慢很多
    while (iter < maxIter):
        # 记录该alpha是否已经进行优化
        alphaPairsChanged = 0
        # 一对x(i),y(i)对应一个alphas(i)
        # 遍历每一个alphas(相当于选择了第一个alphas)，再在剩下的alphas中寻找第二个alphas
        for i in range(m):
            # 根据当前的所有alphas计算该样本xi的预测值yi
            # yi计算公式 y[i] = w^Tx[i]+b
            # w计算公式  w = Σ(1~n) a[n]*label[n]*x[n]
            fXi = float(multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[i, :].T)) + b
            # 计算误差Ei 预测值减真实值
            Ei = fXi - float(labelMat[i])
            # 满足两个优化条件则alphas[i]才能进行优化
            # (1)如果误差很大说明需要优化 : yi预测错误的概率 > 容错率
            #       yi预测错误的期望 : labelMat[i]*Ei
            #       容错率 : toler
            # (2)alphas要在边界之内 : 0<=alphas[i]<=C
            #       因为优化方式是增加一个alphas和降低一个alphas，alphas=0则不能减少，alphas=C则不能再增加
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                # 选第二个alpha
                j = selectJrand(i, m)
                # 计算yj，以及yj预测错误的期望
                fXj = float(multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[j, :].T)) + b
                Ej = fXj - float(labelMat[j])
                # 保存优化之前的alphas(i)和alphas(j)
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()
                # 计算L和H 参考李航统计学方法126 L和H是alphas[j]对角线段 端点的界，必须满足 L <= alphas[j] <= H
                # y(i)和y(j)异侧
                if (labelMat[i] != labelMat[j]):
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                # y(i)和y(j)同侧
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                # 如果相同，就没发优化了，进入到下一个alphas[i]过程中
                if L == H:
                    print("L==H")
                    continue
                # 计算 eta(即alphas[j]的最优调整量) eta = Kii + Kjj - 2Kij 参考李航统计学方法p127
                eta = 2.0 * dataMatrix[i, :]*dataMatrix[j, :].T - dataMatrix[i, :]*dataMatrix[i, :].T - dataMatrix[j, :]*dataMatrix[j, :].T
                if eta >= 0:
                    print("eta>=0")
                    continue
                # 计算调整后的alphas[j]值  参考李航统计学方法p127
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
                # 将新alphas[j]调整到L和H之间
                alphas[j] = clipAlpha(alphas[j], H, L)
                # alpha[j]如果是轻微改变的话，alphas[i]就没必要再调整了
                if (abs(alphas[j] - alphaJold) < 0.00001):
                    print("j not moving enough")
                    continue
                # 调整alphas[i]，alphas[i]和alphas[j]调整量一样只是一个增加一个减小
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
                # 在对alpha[i], alpha[j] 进行优化之后，更新一下b
                # b为位移项决定了超平面与原点的距离
                # w= Σ[1~n] ai*yi*xi
                # y = w'x + b
                # b = yj- Σ[1~n] ai*yi(xi*xj)
                # b1 - b = (y1-y) - Σ[1~n] yi*(a1-a)*(xi*x1)
                b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i, :]*dataMatrix[i, :].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i, :]*dataMatrix[j, :].T
                b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i, :]*dataMatrix[j, :].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j, :]*dataMatrix[j, :].T
                if (0 < alphas[i]) and (C > alphas[i]):
                    b = b1
                elif (0 < alphas[j]) and (C > alphas[j]):
                    b = b2
                else:
                    b = (b1 + b2)/2.0
                alphaPairsChanged += 1
                print("iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged))
        # 在for循环外，检查alpha值是否做了更新，如果在更新则将iter设为0后继续运行程序
        # 知道更新完毕后，iter次循环无变化，才推出循环。
        if (alphaPairsChanged == 0):
            iter += 1
        else:
            iter = 0
        print("iteration number: %d" % iter)
    return b, alphas


'''
基于alphas和xi计算w，算出w就可以算出yi
yi计算公式 y[i] = w^Tx[i]+b
w计算公式  w = Σ(1~n) a[n]*label[n]*x[n]
'''
def calcWs(alphas, dataArr, classLabels):
    X = mat(dataArr)
    labelMat = mat(classLabels).transpose()
    m, n = shape(X)
    w = zeros((n, 1))
    for i in range(m):
        w += multiply(alphas[i] * labelMat[i], X[i, :].T)
    return w


'''
画出分界线和支持向量
'''
def plotfig_SVM(xMat, yMat, ws, b, alphas):
    '''
    xMat:
    yMat:
    ws:  W
    b:   B
    alphas:  alpha
    '''
    xMat = mat(xMat)
    yMat = mat(yMat)
    # b原来是矩阵，先转为数组类型后其数组大小为（1,1），所以后面加[0]，变为(1,)
    b = array(b)[0]
    fig = plt.figure()
    ax = fig.add_subplot(111)
    # 注意flatten的用法
    ax.scatter(xMat[:, 0].flatten().A[0], xMat[:, 1].flatten().A[0])
    # x最大值，最小值根据原数据集dataArr[:, 0]的大小而定
    x = arange(-1.0, 10.0, 0.1)
    # 根据x.w + b = 0 得到，其式子展开为w0.x1 + w1.x2 + b = 0, x2就是y值
    y = (-b-ws[0, 0]*x)/ws[1, 0]
    ax.plot(x, y)

    for i in range(shape(yMat[0, :])[1]):
        if yMat[0, i] > 0:
            ax.plot(xMat[i, 0], xMat[i, 1], 'cx')
        else:
            ax.plot(xMat[i, 0], xMat[i, 1], 'kp')

    # 找到支持向量，并在图中标红
    for i in range(100):
        if alphas[i] > 0.0:
            ax.plot(xMat[i, 0], xMat[i, 1], 'ro')
    plt.show()


if __name__ == "__main__":
    dataArr, labelArr = loadDataSet('svmTestSet.txt')
    # 计算b和alphas
    b, alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40)
    print('/n/n/n')
    print('b=', b)
    print('alphas[alphas>0]=', alphas[alphas > 0])
    print('shape(alphas[alphas > 0])=', shape(alphas[alphas > 0]))
    for i in range(100):
        if alphas[i] > 0:
            print(dataArr[i], labelArr[i])
    # 计算w
    ws = calcWs(alphas, dataArr, labelArr)
    plotfig_SVM(dataArr, labelArr, ws, b, alphas)